home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpklib.zip / SQRTS.FOR < prev    next >
Text File  |  1984-01-06  |  15KB  |  454 lines
  1.       SUBROUTINE SQRTS(LUNIT)
  2. C     LUNIT IS THE OUTPUT UNIT NUMBER
  3. C
  4. C     TESTS
  5. C        SQRDC,SQRSL
  6. C
  7. C     LINPACK. THIS VERSION DATED 08/14/78 .
  8. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  9. C
  10. C     SUBROUTINES AND FUNCTIONS
  11. C
  12. C     LINPACK SQRDC,SQRSL
  13. C     EXTERNAL SQRBM,SQRFM,SQRRX,SXGEN,SMACH
  14. C     FORTRAN AMAX1,ABS,FLOAT
  15. C
  16. C     INTERNAL VARIABLES
  17. C
  18.       INTEGER LUNIT
  19.       INTEGER CASE,LDX,N,P,PD2,PD2P1,NCASE,JPVT(25)
  20.       INTEGER I,INFO,J,JJ,JJJ,L
  21.       REAL BESTAT,RMAX,RSTAT,XBMAX,XBSTAT,ERRLVL,BSTAT,FSTAT
  22.       REAL BETA(25),QRAUX(25),QY(25),RSD(25),S(25),WORK(25)
  23.       REAL X(25,25),XBETA(25),XX(25,25),Y(25),Y1(25),Y2(25)
  24.       REAL SMACH,TINY,HUGE
  25.       LOGICAL NOTWRT
  26.       LDX = 25
  27.       NCASE = 9
  28.       ERRLVL = 100.0E0
  29.       NOTWRT = .TRUE.
  30.       TINY = SMACH(2)
  31.       HUGE = SMACH(3)
  32.       WRITE (LUNIT,600)
  33.       DO 550 CASE = NCASE, NCASE
  34.          WRITE (LUNIT,560) CASE
  35.          XBSTAT = 0.0E0
  36.          BESTAT = 0.0E0
  37.          GO TO (10, 100, 170, 240, 300, 330, 380, 430, 470), CASE
  38. C
  39. C        WELL CONDITIONED LEAST SQUARES PROBLEM.
  40. C
  41.    10    CONTINUE
  42.             WRITE (LUNIT,640)
  43.             N = 10
  44.             P = 4
  45. C
  46. C           GENERATE A MATRIX X WITH SINGULAR VALUES 1. AND .5.
  47. C
  48.             PD2 = P/2
  49.             PD2P1 = PD2 + 1
  50.             DO 20 L = 1, PD2
  51.                S(L) = 1.0E0
  52.    20       CONTINUE
  53.             DO 30 L = PD2P1, P
  54.                S(L) = 0.5E0
  55.    30       CONTINUE
  56.             CALL SXGEN(X,LDX,N,P,S)
  57.             IF (NOTWRT) GO TO 40
  58.                WRITE (LUNIT,610)
  59.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  60.    40       CONTINUE
  61. C           THE OBSERVATION VECTOR IS Y = Y1 + Y2, WHERE Y1 IS
  62. C           THE VECTOR OF ROW SUMS OF X AND Y2 IS ORTHOGONAL TO
  63. C           THE COLUMNS OF X.
  64. C
  65.             DO 60 I = 1, N
  66.                Y1(I) = 0.0E0
  67.                Y2(I) = 2.0E0*TINY
  68.                IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)*TINY
  69.                DO 50 J = 1, P
  70.                   X(I,J) = X(I,J)*TINY
  71.                   Y1(I) = Y1(I) + X(I,J)
  72.                   XX(I,J) = X(I,J)
  73.    50          CONTINUE
  74.                Y(I) = Y1(I) + Y2(I)
  75.    60       CONTINUE
  76. C
  77. C           REDUCE THE MATRIX.
  78. C
  79.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  80.             IF (NOTWRT) GO TO 70
  81.                WRITE (LUNIT,610)
  82.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  83.                WRITE (LUNIT,620)
  84.                CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
  85.    70       CONTINUE
  86. C
  87. C           COMPUTE THE STATISTICS FOR THE FOWARD AND BACK
  88. C           MULTIPLICATIONS.
  89. C
  90.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  91.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  92. C
  93. C           SOLVE THE LEAST SQUARES PROBLEM.
  94. C
  95.             CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
  96.      *                 INFO)
  97. C
  98. C           COMPUTE THE LEAST SQUARES TEST STATISTICS.
  99. C
  100.             RSTAT = 0.0E0
  101.             RMAX = 0.0E0
  102.             XBSTAT = 0.0E0
  103.             XBMAX = 0.0E0
  104.             DO 80 I = 1, N
  105.                RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
  106.                RMAX = AMAX1(RMAX,ABS(Y2(I)))
  107.                XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
  108.                XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
  109.    80       CONTINUE
  110.             RSTAT = RSTAT/(SMACH(1)*RMAX)
  111.             XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
  112.             BESTAT = 0.0E0
  113.             DO 90 J = 1, P
  114.                BESTAT = AMAX1(BESTAT,ABS(1.0E0-BETA(J)))
  115.    90       CONTINUE
  116.             BESTAT = BESTAT/SMACH(1)
  117.             WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
  118.          GO TO 540
  119. C
  120. C        4 X 10 MATRIX.
  121. C
  122.   100    CONTINUE
  123.             WRITE (LUNIT,650)
  124.             N = 4
  125.             P = 10
  126.             PD2 = P/2
  127.             PD2P1 = PD2 + 1
  128.             DO 110 L = 1, PD2
  129.                S(L) = 1.0E0
  130.   110       CONTINUE
  131.             DO 120 L = PD2P1, P
  132.                S(L) = 0.5E0
  133.   120       CONTINUE
  134.             CALL SXGEN(X,LDX,N,P,S)
  135.             IF (NOTWRT) GO TO 130
  136.                WRITE (LUNIT,610)
  137.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  138.   130       CONTINUE
  139.             DO 150 I = 1, N
  140.                DO 140 J = 1, P
  141.                   XX(I,J) = X(I,J)
  142.   140          CONTINUE
  143.   150       CONTINUE
  144.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  145.             IF (NOTWRT) GO TO 160
  146.                WRITE (LUNIT,610)
  147.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  148.                WRITE (LUNIT,620)
  149.                CALL SARRAY(QRAUX,N,N,1,-4,LUNIT)
  150.   160       CONTINUE
  151.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  152.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  153.             WRITE (LUNIT,570) FSTAT,BSTAT
  154.          GO TO 540
  155. C
  156. C        PIVOTING AND OVERFLOW TEST.  COLUMNS 1, 4, 7 FROZEN.
  157. C
  158.   170    CONTINUE
  159.             WRITE (LUNIT,670)
  160.             N = 10
  161.             P = 3
  162.             S(1) = 1.0E0
  163.             S(2) = 0.5E0
  164.             S(3) = 0.25E0
  165.             CALL SXGEN(X,LDX,N,P,S)
  166.             P = 9
  167.             DO 190 I = 1, N
  168.                DO 180 JJJ = 1, 3
  169.                   J = 4 - JJJ
  170.                   JJ = 3*J
  171.                   X(I,JJ) = HUGE*X(I,J)
  172.                   X(I,JJ-1) = X(I,JJ)/2.0E0
  173.                   X(I,JJ-2) = X(I,JJ-1)/2.0E0
  174.   180          CONTINUE
  175.   190       CONTINUE
  176.             IF (NOTWRT) GO TO 200
  177.                WRITE (LUNIT,610)
  178.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  179.   200       CONTINUE
  180.             DO 220 J = 1, P
  181.                JPVT(J) = 0
  182.                DO 210 I = 1, N
  183.                   XX(I,J) = X(I,J)
  184.   210          CONTINUE
  185.   220       CONTINUE
  186.             JPVT(1) = -1
  187.             JPVT(4) = -1
  188.             JPVT(7) = -1
  189.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  190.             IF (NOTWRT) GO TO 230
  191.                WRITE (LUNIT,610)
  192.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  193.                WRITE (LUNIT,620)
  194.                CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  195.   230       CONTINUE
  196.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  197.             CALL SQRRX(XX,LDX,N,P,JPVT)
  198.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  199.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  200.             WRITE (LUNIT,570) FSTAT,BSTAT
  201.          GO TO 540
  202. C
  203. C        25 BY 25 MATRIX
  204. C
  205.   240    CONTINUE
  206.             WRITE (LUNIT,680)
  207.             N = 25
  208.             P = 25
  209.             DO 250 I = 1, 25
  210.                S(I) = 2.0E0**(1 - I)
  211.   250       CONTINUE
  212.             CALL SXGEN(X,LDX,N,P,S)
  213.             IF (NOTWRT) GO TO 260
  214.                WRITE (LUNIT,610)
  215.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  216.   260       CONTINUE
  217.             DO 280 J = 1, P
  218.                JPVT(J) = 0
  219.                DO 270 I = 1, N
  220.                   XX(I,J) = X(I,J)
  221.   270          CONTINUE
  222.   280       CONTINUE
  223.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  224.             IF (NOTWRT) GO TO 290
  225.                WRITE (LUNIT,610)
  226.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  227.                WRITE (LUNIT,620)
  228.                CALL SARRAY(QRAUX,P,P,1,-10,LUNIT)
  229.   290       CONTINUE
  230.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  231.             CALL SQRRX(XX,LDX,N,P,JPVT)
  232.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  233.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  234.             WRITE (LUNIT,570) FSTAT,BSTAT
  235.          GO TO 540
  236. C
  237. C        MONOELEMENTAL MATRIX.
  238. C
  239.   300    CONTINUE
  240.             WRITE (LUNIT,690)
  241.             N = 1
  242.             P = 1
  243.             X(1,1) = 1.0E0
  244.             IF (NOTWRT) GO TO 310
  245.                WRITE (LUNIT,610)
  246.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  247.   310       CONTINUE
  248.             XX(1,1) = 1.0E0
  249.             JPVT(1) = 0
  250.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  251.             IF (NOTWRT) GO TO 320
  252.                WRITE (LUNIT,610)
  253.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  254.                WRITE (LUNIT,620)
  255.                CALL SARRAY(QRAUX,P,P,1,-1,LUNIT)
  256.   320       CONTINUE
  257.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  258.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  259.             WRITE (LUNIT,570) FSTAT,BSTAT
  260.